bellman.m

Calculates LIFO equilibrium continuation values using Bellman iteration and plots the results.

Contents

Set up the problem's parameters.

This section first checks to see if the parameters of the problem are defined in the structure bellmanARG . If this does not exist, then it defines them.

if exist('bellmanARG','var')

    %Parameters for constructing the Markov chain.
    markovARG=bellmanARG.markovARG;

    %Producer's surplus per customer as a function of N.
    piF=bellmanARG.piF;

    %Sunk entry cost and per-period fixed cost.
    phi=bellmanARG.phi;
    kappa=bellmanARG.kappa;

    %Discount factor
    beta=bellmanARG.beta;

else

    %Parameters for approximating the innovation distribution
    approximateARG.F=@(x) (1+erf(x/sqrt(2)))./2;
    approximateARG.Frange=3.0;
    approximateARG.k=51;

    %Parameters for constructing the Markov chain.
    markovARG.approximateARG=approximateARG;
    markovARG.rho=1.0;
    markovARG.sigma=0.30;

    markovARG.omegaCenter=0.0;
    markovARG.omegaStep=0.01;
    markovARG.omegaWidth=3;

    %Parameters describing profits and
    %discounting.

    k=4;
    piF = @(x) k*ones(size(x));     %Producer surplus per customer
    kappa=1.75;                     %per-period fixed cost
    beta=1.05^(-1);                 %5 percent annual interest rate.
	phi = @(N) 0.25*beta/(1-beta);  %Sunk cost of entry

end

Create Markov chain if required.

The results are stored in Pi and omega.

if ~exist('Pi','var')
    markov
end

Calculate upper bound on firm numbers.

nCheck=find((max(exp(omega))*piF(1:1:100)./(1:1:100)>=kappa),1,'last');

Solve individual entrants' problems

% Store the individual firms' value functions in a three-dimensional array.
% The final dimension gives the rank of the firm.
% The second dimension gives the number of firms presently active.
% The first dimension indexes the demand states.

vFuncs=zeros(nCstates,nCheck,nCheck);


for R=nCheck:-1:1

    %Initialize the value function.
    v=zeros(nCstates,nCheck-R+1);

    %Initialize the transition rule for N if necessary.
    if R==nCheck;
        nPrime=nCheck*ones(nCstates,1);
    end

    %Create nPrimeIndex, which determines which column of v(:,:) contains the relevant continuation value.
    if R==nCheck
        nPrimeIndex=ones(nCstates,nCheck-R+1);
    else
        nPrimeIndex=zeros(nCstates,nCheck-R+1);
        for i=R:1:nCheck;
            nPrimeIndex=nPrimeIndex+(nPrime==i)*(1+(i-R));
        end
    end

    % Bellman Iteration
    tol=1.0e-7;
    vDistance=1;

    while vDistance>=tol

        %Use nprime to assign the appropriate continuation value.
        Piv=Pi*v;
        Ev=zeros(nCstates,nCheck-R+1);
        for i=1:nCstates;
            Ev(i,:)=Piv(i,nPrimeIndex(i,:));
        end

        vPrime=beta*(Pi*exp(omega')*ones(1,nCheck-R+1)).*piF(nPrime)./nPrime-beta*kappa+beta*Ev;

        %Take the maximum of the expected continuation value and zero.
        vPrime=max(vPrime,0);

        %Measure the distance between the initial and new value functions.
        vDistance=max(max(abs(vPrime-v)));

        %Update the value function.
        v=vPrime;

    end

    %Store the value function
    if R>1; %In this case, we need to pad the matrix of values with NaNs
        vv=[NaN(nCstates,R-1) v];
    else
        vv=v;
    end

    vFuncs(:,:,R)=vv;

    %Update nprime for the problem of a firm with rank R-1

    %Take away one firm if exit is optimal.
    nPrime=nPrime-(v==0);

    %Add a column to nPrime for the initial state N=R-1 (The Matlab editor complains about the speed consequences of
    %this assigment. For a big problem where the speed drag from copying nPrime to a new location in memory is too much
    %to pay, preallocate nPrime to a nCstates x nCheck+1 array and only use the relevant column in the Bellman equation
    %above.)
    nPrime=[(R-1)*ones(nCstates,1) nPrime];

    %Add firms to the new column in the states where entry is optimal.
    for i=R:1:nCheck
        nPrime(:,1)=nPrime(:,1)+(squeeze(vFuncs(:,i,i))>phi(i));
    end

end

Create final version of nPrimeIndex.

nPrimeIndex=zeros(nCstates,nCheck+1);
for i=1:1:nCheck;
    nPrimeIndex=nPrimeIndex+(nPrime==i)*(1+i);
end

Determine the ergodic set for N.

minN=min(min(nPrime));
maxN=max((max(nPrime).*(max(nPrime)>(0:1:nCheck)))); %Highest value of N transited to from a lower value.

%Trim the matrix nPrime to eliminate states outside of the ergodic set.
nPrime=nPrime(:,minN+1:maxN+1);
nPrimeIndex=nPrimeIndex(:,minN+1:maxN+1)-minN+1;

%Trim the value functions to eliminate values
%of N above the ergodic set.
vFuncs=vFuncs(:,1:maxN,1:maxN);

Plot the value functions.

This section creates a raw inelegant plot which can be used to verify that the value functions satisfy the restrictions of theory.

if ~exist('bellmanARG','var')

    figure(1)
    set(gcf,'Units','inches','Position',[0 0 4 6]);
    for i=1:maxN;
        subplot(maxN,1,i);
        vi=squeeze(vFuncs(:,:,i));
        plot(omega'*ones(1,maxN-i+1),vi(:,i:maxN)); %This line plots the value function's many ``branches'' all at one go.
        hold
        plot(omega',phi(i)*ones(size(omega)));      %This line adds a horizontal line to mark the entry cost.
    end
end
Current plot held
Current plot held
Current plot held
Current plot held